function bisect(f, a, b, tol)
  implicit none
  real :: bisect
  real, external :: f
  real, intent(in) :: a, b, tol
  real :: c, fa, fb, fc, aa, bb
  integer :: k
!  write(*,*) 'hfhhlsfls'
  fa = f(a)
  fb = f(b)
  aa = a
  bb = b
  write(*,*) fa, fb
  if(fa*fb >= 0) then
     write(*,*) 'fa*fb < 0 not satisfied'
  end if

  do while((bb-aa)/2 > tol)
     c = (aa+bb)/2
 !    write(*,*) c
     fc = f(c)
     if(fc == 0.0) exit
     if(fc*fa < 0.0) then
        bb = c
        fb = fc
     else
        aa = c
        fa = fc
     end if
  end do

  bisect = (aa+bb)/2
end function bisect


function f(x)
  implicit none
  real :: f
  real, intent(in) :: x
  f = x**3 + x -1.0
end function f


program test_bisect
  implicit none
  real, external :: f
  real :: bisect
  write(*,*) bisect(f, 0.0, 1.0, 0.00005)
end program test_bisect
